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Abstract. The wavefronts associated with a one-dimensional combustion model with Arrhe- 
nius kinetics and no heat loss are analyzed within the high Lewis number perturbative limit. This 
situation, in which fuel diffusivity is small in comparison to that of heat, is appropriate for highly 
dense fluids. A formula for the wavespeed is established by a non-standard application of Melnikov's 
method and slow manifold theory from dynamical systems, and compared to numerical results. A 
simple characterization of the wavespeed correction is obtained: it is proportional to the ratio be- 
tween the exothermicity parameter and the Lewis number. The perturbation method developed 
herein is also applicable to more general coupled reaction-diffusion equations with strongly differing 
diffusivities. The stability of the wavefronts is also tested using a numerical Evans function method. 
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1. Introduction. In this article, we study the wavespeed of a combustion wave- 
front along a one-dimensional medium. This is a fundamental idealized problem 
towards understanding how flame fronts propagate, and therefore has received a con- 
siderable amount of attention. There are several (non-dimensional) parameters of 
importance: the Lewis number Le, the exothermicity parameter /?, and the heat loss 
parameter £. The first of these, the Lewis number, measures the relative importance 
of fuel diffusivity in comparison to that of heat. The exothermicity (3 is the ratio of 
the activation energy to the heat of reaction. The structure of the governing equa- 
tions is such that an infinite Lewis number is considerably easier to deal with than 
allowing for fuel diffusivity. Many studies of this "solid" regime appear in the lit- 
erature [5,7,28,36,37], and also the "gaseous" regime Le » 1 has been frequently 
studied because of a symmetry in the equations [7,20,24,37,39]. Usually, the heat 
loss is neglected in these "adiabatic" studies. In several of these articles [28,36,37] 
the condition /3 ^> 1 is essential to the wavespeed and stability analysis. The case 
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j3 <C 1 has also been studied [8], in which a perturbative method is used to model 
the temperature. The bifurcation structure with respect to the heat loss parameter 
I is addressed in [34], which obtains a stability diagram with respect to t and the 
wavespeed. 

We note that the limit of small fuel diffusivity (large, but not infinite, Lewis 
number) has not received much attention, perhaps because of the singularity of this 
limit in the governing equations. Yet this limit may be argued to be particularly 
appropriate for very high density fluids burning at high temperatures, such as would 
occur, for example, in the burning of toxic wastes at supercritical temperatures [25]. 
Even for solids, some mass diffusivity is to be expected at very high temperatures, 
particularly in the reaction zone in which liquification may occur. In [27], the mass 
diffusivity is modeled by an Arrhenius temperature dependence, which would result 
in a large effective Lewis number in certain situations (such as when the [scaled] 
adiabatic flame temperature is small in comparison to the activation energy for mass 
diffusion). It is this very large Lewis number limit which we study in this article, 
without restricting f3. We do a detailed analysis of the wavespeed of combustion 
waves which can be supported. We also verify the linear stability of such wavefronts 
using an Evans function technique. 

The model we use is for a premixed fuel in one dimension, with no heat loss 
and with an Arrhenius law for the reaction rate. These combustion dynamics can be 
represented in non-dimensional form by [5,8,20,24,28,34,36,37,39] 

du d 2 u 



(1.1) 



at = d^ + ye ~ 1/u 



dt Le dx 2 

Here, u{x,t) is the temperature, and y(x,t) the fuel concentration, at a point x at 
time t. The parameters (3 and Le are as described earlier. We are neglecting heat loss 
(had we included it, an additional term — t (u — u a ) for some ambient temperature 
u a would be necessary on the right-hand side of the it equation in {HI}). This one- 
dimensional model is also applicable to combustion in cylinders [30], with u and y 
being cross-sectionally averaged quantities in this case. See also [6,7,19,26] for closely 
related governing equations. The non-dimensionalization leading to (|1.1[) ensures that 
the cold boundary problem is circumvented (see [36] for a discussion). Since the Lewis 
number will be assumed large, set e = 1/Le with 0<e<l. This small e limit clearly 
constitutes a singular perturbation in (|1.1[) . 

This article analyzes (jl.lj) as follows. In Section [21 we determine the wavespeed 
as a function of f3 and e. We initially consider the situation where Le = oo (Sec- 
tion [2T}, since this wavespeed is relevant to our subsequent perturbative analysis for 
1 <C Le < oo (Sections 12.21 12.31 and 12. 4| ). While the infinite Lewis number situation is 
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well-studied, we are able to empirically determine a simple exponential formula for the 
wavespeed as a function of (3. The case 1 <C Le < oo is initially examined numerically 
in Section 12. 2[ in which we obtain a method for computing the wavespeed. In the 
subsequent sections, we establish a theoretical estimate for the wavespeed with the 
help of two suitably modified tools from dynamical systems theory: a slow manifold 
reduction, and Melnikov's method. In Section 12. 3[ we reduce the dimensionality of 
the problem using a slow manifold reduction argument. This enables us in Section HOI 
to utilize a non-standard adaptation of Melnikov's method to find a theoretical es- 
timate for the wavespeed. (This new technique is adaptable to other situations in 
which the wavespeed correction due to the presence of a small parameter is needed.) 
Our asymptotics enable the determination of a remarkably simple formula for the 
wavespeed, which is accurate for all fj values (and not restricted to the "usual" large 
(3 limit). Essentially, we find that the relative wavespeed correction in going from 
infinite to large Lewis number is proportional to (/3/Le). 

A brief stability analysis of the wavefronts is given in Section RJ] Having described 
the Evans function approach to stability in Section [3l we compute the Evans function 
for high Lewis number combustion wavefronts using an exterior algebra [2, 16, 23, 
40]. We note that in [20], an exterior algebra method has been successfully used 
to numerically investigate stability of wavefronts in combustion systems. A detailed 
stability analysis in the (3 - Le -1 plane is given therein for the system (|1.1[) . As with 
infinite Lewis number fronts (see [5,14,28,37]), [20] show that stability occurs for small 
(3, but that as (3 is increased, a Hopf bifurcation leads to an oscillatory instability. The 
(3 and Le values we test give results consistent with the stability boundary determined 
in [20]. Thus, stability properties remain essentially unaltered despite the singularity 
in the limit Le — > oo. 

2. Wavespeed analysis. We seek wavefronts which travel in time, and hence 
set u(x,t) = it(£) and y(x,t) = j/(£), where £ = x — ct and c is the traveling wave 
speed. Under this ansatz, (|1.1|) reduces to 



2.1. Wavefront for Le = oo. Set e = in (|2.1| . Upon defining the new variable 
v = u' , the dynamics can be represented by a three-dimensional first-order system 



(2.1) 












v! = V 



(2.2) 



< 



V 



c v 



ye 



l/u 
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Fig. 2.1. Projection onto the (u,v)-plane of trajectories of i2.2\) lying on different planes 
H c = c. Here, (3 = 1, and the three curves correspond to c = 0.5 (dotted), 0.5707 (solid) and 0.7 
( dashed). 



The system (|2.2[) possesses a conserved quantity 

(2.3) H c (u,v,y) = (3v + (3cu + cy , 

since it is verifiable that dH c /d^ = along trajectories of (|2.2|) . Thus, motion is 
confined to planes defined by H c (u, v, y) = constant. Now, for a wavefront, we require 
that (u, v, y) — > (0, 0, 1) as £ — * oo; this corresponds to the region in which fuel is not 
yet burnt (and remains at its maximum non-dimensional concentration of one) and 
the temperature (and its variation) is still zero. This point lies on H c (u,v,y) = c, 
giving a well-known conservation relation [37]. At the other limit £ — > —oo, the fuel is 
completely burnt, and has reached a steady temperature, and so (it, v,y) — * (ub, 0, 0), 
where the temperature ub is to be determined. Utilizing H c (ub,Q,0) — c, we find 
that ub = is necessary; see also [8,20,37] for alternative ways to obtain this value. 

Thus, we seek a heteroclinic solution of (|2.2p . which progresses between the fixed 
points (1//3, 0, 0) and (0, 0, 1), and is confined to the plane v + (3 cu + cy = c; i.e., 
the fuel concentration obeys 

(2.4) y=l-Pu--v 
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Fig. 2.2. Variation of the wavespeed c with j3: open circles (numerical results); unbroken curve 
(empirical curve 12. 61) ); dotted curve (exp(— 0.5/3), as obtained in [7, 28, 36]). 



at all values of £. Considering (|2.2[) under this restriction, we obtain 

vf = v 

v' = -cv- (l-(3u-^v) e- 1/u 

This is effectively a projection of the flow on the particular invariant plane H c (u, v, y) 
= c onto the (it, w)-plane. Any value of c for which a heteroclinic connection exists 
between (1/(3, 0) and (0, 0) is a permitted speed for the wavefront. 

The unstable eigen-direction of the point (1/(3, 0) is (— c, —(3e~P), and we deter- 
mine c numerically by shooting along this direction, and attempting to match up 
with a trajectory approaching the origin. In Figure 12. 1 1 we show several numerically 
computed trajectories of this form, for different values of c, where we have chosen 
(3 = 1. Note that this is not a standard (it, u)-phase space for (|2.5jl . since each curve 
corresponds to a different value of the parameter c. Rather, it is a projection onto 
the (it, i>)-plane of specialized trajectories from the invariant planes H c (u,v,y) = c 
of the three-dimensional system (|2.2p . The one trajectory which makes the required 
connection lies in the invariant plane corresponding to c = 0.5707. The determina- 
tion of this c value was obtained by making incremental adjustments of c until an 
appropriate connection is obtained. 
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We use this technique to numerically compute the wavespeeds for various values 
of the fuel parameter [3, and illustrate this dependence in Figure 12.21 The wavespeed 
decays with (3. For fuels with larger (3 (poor fuels) , the energy resulting from the reac- 
tion is insufficient to quickly activate combustion in nearby material, and combustion 
fronts propagate slowly. The data fits the exponential curve 

(2.6) c(/3) = 0.927 e- a486/3 

with correlation p > 0.9999. Equation (|2.6p therefore provides an empirically deter- 
mined formula of excellent accuracy, for the speed of a wavefront in perfectly solid 
adiabatic one-dimensional media. This result is close (and consistent with) a variety 
of sources: exp(— 0.5/3) is quoted in [36] for the small (3 limit; this same value is given 
as an upper bound in [7], and also implied in eq. (10) in [28] using a large (3 limit 
within a discontinuous front approximation. See Figure [2~2l for a comparison with our 
results. 

The structure of the temperature front is illustrated in Figure |2~31 for (3—1 (solid 
curve, left scale) and (3 = 3 (dashed curve, right scale), demonstrating that larger 
(3 fronts have a broader preheat layer preceding the front. Note that the preheat 
zone differs from the reaction zone [24]. The latter shrinks with increasing {3 [13]. 
Specifically, the reaction zone as a fraction of the preheat zone is 0(1/ (3) [24,28]. The 
reaction zone is well localized near the region of greatest temperature change [24], 
and is not immediately identifiable in temperature profiles as in Figure 12.31 Indeed, 
the increase in size of the preheat layer with (3 supports the 0(1/ [3) expectation for 
the ratio between the reaction and the preheat zones. 

2.2. Wavespeed for 1 <C Le < oo. When the Lewis number is not infinite, but 
large, e is small, and weak fuel diffusion needs to be permitted. This is a singular hmit 
in (jl.ip and (|2.1[) . and as a consequence has been hardly examined in the literature. 
By defining v = v! as before, but now also z = y', the governing equations (|2.1j) can 
be represented as a four-dimensional system 

vl = v 



(2.7) 



v' = -cv-ye x l u 
y' = z 

z' = - (-cz + f3ye~ 1/j 
e V 

This is reducible to three-dimensions: the quantity 

G e c (u,v,y, z) = [3v + f3cu + cy + ez 

can be verified to be a conserved quantity of (|2.7| . Hence, flow is confined to the 
invariant three-dimensional surfaces G% = constant. Now, we seek a wavefront solution 
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Fig. 2.3. Temperature front at f3 = 1 (solid, left scale) and [3 = 3 (dashed, right scale) 



which goes from (u,v,y,z) = (ub, 0,0,0) to a value (0,0,1,0), and we find that 
GJ(u,w,y,2) = c, and ug — 1/(3 as before. The three-dimensional invariant surface 
on which both points lie is 

z = — (c — (3 v — j3 cu — cy) . 
e 

The dynamics of (|2.7p on this surface can be projected to the three variables (u, v, y), 
such that 



(2.8) 



= —cv — ye 



-i/i 



- (c — (3 v — (3 cu — cy) 



We seek the value of c which permits a heteroclinic connection from (u,v,y) = 
(1//3, 0, 0) to (0,0,1). The former point (corresponding to £ = — oo) has only one 
positive eigenvalue, given by (-c+ v/c 2 +Ae(3e-< 3 )/{2e). For small e, we "shoot" in 
the eigen-direction corresponding to this, with an initial guess of the wavespeed given 
by (|2.6|) . Thereafter, as in the previous section, we adjust c until we obtain a solution 
which approaches the point (0, 0, 1) as £ — -> oo. We do this numerically by considering 
the conditions cy + ez = c and v + cu = which the front must obey at £ = +oo, 



8 BALASURIYA, GOTTWALD, HORNIBROOK AND LAFORTUNE 
0.58 i 1 1 1 1 1 




0.52 I ' ' ' ' ' ' 1 

0.05 0.1 0.15 0.2 0.25 0.3 

e 

Fig. 2.4. Numerically computed wavespeed variation with e for f3 = 1 (the crosses). The dashed 
line is the theoretical approximation for (3 = 1 obtained from the methods of Section \2.4\ 

and using a root-finding algorithm to adjust c. For a fixed value (3 = 1, we illustrate 
how the wavespeed c varies with e in Figure [2^41 with the crosses. The dashed curve 
in Figure 12^41 is a analytical/numerical approximation we obtain for the wavespeed in 
terms of an explicit formula (|2.17p . The next two sections describe how we obtain 
this formula. 

We notice that c decreases as we increase e, that is, when we decrease the Lewis 
number. Now, in dimensional form Le = K,/(pc p D), where p,n,c p and D are the 
respectively the density, thermal conductivity, specific heat capacity and molecular 
diffusivity of the fuel [6,8,30,34,37]. Increasing e is equivalent to increasing the relative 
importance of D, p and c p in relation to k. Reducing k obviously decreases the ability 
of heat to move, and hence the combustion speed. Higher densities result in increased 
fuel mass in each location, which means more heat is needed in a given area to ignite 
all of the fuel before the wave moves on. Fuels with increased c p require more heat to 
increase the temperature by the a specified amount. Finally, increasing D increases 
the transport of burnt fuel into the unburnt region and vice-versa, interfering with 
front propagation. 

We computed the changes to the wavefront profile (akin to Figure 12. 3|) when e 
is changed (not shown). We verified the obvious physical conclusion that the fuel 
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Fig. 2.5. The hyperbolic invariant manifold (a) So for \Z. 10p . and (b) <S £ for (2. 91) 

concentration front becomes less steep when e is increased from zero. 

2.3. Slow manifold reduction. We now show that in the limit of small e, it 
is possible to further reduce the system (|2.8p to a fiuo-dimensional flow on a slow 
manifold. We begin with l|2.8[) . and note that there are two "time" -scales in this 
singularly perturbed system, where we use "time" loosely to mean the independent 
variable £. We therefore adopt the standard dynamical systems trick of defining a 
new independent variable 77 = £/e to elucidate motion in the fast "time" 77. With a 
dot denoting the rate of change with respect to n, (|2.8p becomes 



in which it is clear that the plane <So defined by c — v — (3 cu — cy = consists 
entirely of fixed points. This is the same plane as defined through H c (u, v.y) = c for 
equation (|2.2[) . on which the interesting behavior occurred for perfectly solid fuels. 
Each fixed point has a one-dimensional stable manifold (in the y-direction), and a 
two-dimensional center manifold, which is Sq. Thus the plane So is invariant and 
normally hyperbolic with respect to (|2. 10|1 : there is exponential contraction towards 
it as illustrated in Figure [2~5T aK 



(2.9) 




In the e = limit, the system collapses to 



(2.10) 
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Upon switching on e and considering the dynamics (|2.9|) , iS perturbs to an invari- 
ant curved entity S e , which is order e away from 1S0 . This is because of the structural 
stability of normally hyperbolic sets [18], which also implies that normal hyperbolicity 
is preserved for small e. Therefore, there is exponential decay of trajectories towards 
iS e on "time"-scales of order 77, as shown in Figure f2~5T b) . Motion on S e occurs at a 
slower rate (on "time" -scales of order £), and hence it is termed a 'slow manifold'. The 
heteroclinic connection we seek lies on <S e , from (u, v, y) = (1/(3, 0, 0) to (0, 0, 1). Since 
S e is invariant, two-dimensional and not parallel to the y-axis, it therefore makes sense 
to project the motion to the (u, w)-plane in order to describe behavior. To elucidate 
this motion, we need to once again return to the original "time" -scale £ - the slow 
time associated with motion on the slow manifold. 

Return to the relationship («(£), v{£), z (£)) — c, which upon differentia- 
tion yields 

/3 v' + (3 c v! + cy' + e z' = , 

and since v! = v and y' — z, 

& 1 a e / 

z = v — pv Z ■ 

c c 

Substituting back into G^,(u, v, y, z) — c, we obtain 

(3v + ficu + cy + e \--v' -f3v + 0(e) ] =c, 



and thus 

1 t'-Bti + e- D + e — 



y = l- -v-/3u + e-^v' + e-v + 0(e z ) . 
c c l c 



Substitution into the v' equation in (|2.7[) or (|2.8[) gives 



v 



Therefore 



'il + e^e-vA = - C v-(l-?-v-(3u + e?-v + 0(e 2 ))e-V u . 



c 2 I \ c c 



v'= fl-e4e- 1 /" 



-rr- i I -V j„ . , l r )< 1 " 



Retaining only 0(e) terms, we obtain the (it, w)-projected approximate equations on 
the slow manifold 



(2 ' n) ' v' = -cv~(l--v-f3u)e- 1/u + e^(l--v-f3u)e- 2 / 1 - 



We will now show how to use these approximate dynamics to predict the correction 
to the wavespeed resulting from the inclusion of the finiteness of the Lewis number. 
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Fig. 2.6. Manifold structure for the Melnikov approach: (a) e = 0, (b) e ^ 

2.4. Perturbative formula for wavespeed. Here, we derive and numerically 
study a formula for the wavespeed correction in going from Le = oo to finite Lewis 
number. Let 

(2.12) c 08, e) = c {(3) + e a (/?) + (e 2 ) , 

where Co is the wavespeed associated with the infinite Lewis number (e = 0) com- 
bustion wavefront. In the spirit of perturbation analysis, we obtain a formula for the 
correction c\((3) purely in terms of the unperturbed (e = 0) wave, using a nontradi- 
tional application of "Melnikov's method" [29] from dynamical systems theory. 

Mclnikov's method is applied most commonly to area-preserving two-dimensional 
systems under time-periodic perturbations [4,21,38]. (Here, once again, £ represents 
"time".) Our system (I2.1ip is not area-preserving, and has a perturbation which is 
independent of the temporal variable. Under these conditions, we describe the method 
applied to the system 

(2.13) z'=f(z)+eg(z) . 

When e = 0, suppose this system possesses a heteroclinic connection between the 
two saddle fixed points a and b as shown in Figure [2T6l a). A heteroclinic connection 
of this sort occurs when a branch of the one-dimensional unstable manifold of a 
coincides with a branch of the stable manifold of b. This heteroclinic trajectory can 
be represented as a solution z = z(£) to (|2.13|) with e = 0. 

Now, for small e > in (|2.13p . the fixed points a and b perturb by 0(e), and re- 
tain their stable and unstable manifolds [18]. However, these need no longer coincide. 
Figure [2~6f b) shows how they can split apart, with the dashed curve showing the orig- 
inal manifold. Let g?(£, e) be a distance measure between these manifolds, measured 
along a perpendicular to the unperturbed heteroclinic drawn at z(— £). The variable £ 
can thus be used to identify the position along the heteroclinic curve (cf. "heteroclinic 
coordinates" of Section 4.5 in [38]). Since d(£,0) = for all £, this distance is Taylor 
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expandable in e in the form 



M(0 



|f(*(-0)l 



0(e 2 



where the scaling factor |/(z(— £))| in the denominator represents the unperturbed 
trajectory's speed at the location f. The quantity M(£) is the "Melnikov function", 
for which an expression is 



(2.14) M(0 



/oo r /-m 
exp - / V • f(z(s))ds 
-OO «/ — £ 



f(4(/i)) Ag(z(/x))d/i. 



where the wedge product between two vectors is defined by (ai,a 2 ) T A (6i,&2) T = 
0162 — a2&i. Obtaining the version (|2.14p requires two adjustments to the standard 
Melnikov approaches [4,21,38]: incorporating the non area-preserving nature of the 
unperturbed flow of (|2. 13|1 . and representing the distance in terms of heteroclinic 
coordinates. Details are provided in Appendix |A"1 We need to ensure the persistence 
of a heteroclinic trajectory in (|2.1ip for e > 0, and thus require d(£, e) = for all £. 
For this to happen for all small e, we therefore need to set M(£) = 0. 

To apply this technique to our system, we begin by writing (|2.11[) in the form 
(|2.13|) . Using the expansion (|2.12p . and utilizing binomial expansions for l/(co + eci), 
we get 



(2.15) 



f3c\e~ 



' = -c v-e 1/u T uv +e\ - < <■■ 
where higher-order terms in e have been discarded, and 

T — 1 



1/u 



pu V . 

CO 



By appropriately identifying f and g from (|2.15[) through comparison with (|2.13p . we 
see that 



(f A g) (u, v) = v ( — c\v 



-i/S 



0e 



-2/u 



-T 

-L 1. 



and V • f = -co + /3e _1 / u /co. Substituting into the Melnikov formula (|2.14j) . and 
setting it equal to zero, we obtain 



exp 



co 



I 

CO 



e" 1 '^ ds 



-C\V 



2 uv 



'0 



d/i = , 



where each of u(/j,) and v(/j,) is evaluated along the e = combustion wave. Notice, 
however, that for this infinite Lewis number combustion wave, (12.41) tells us that the 
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Fig. 2.7. TVie perturbing wavespeed as a function of (3 



fuel concentration = T ulI (/x) for all /i. Therefore 



(2.16) Cl (/3)=/3 



exp 



Co 



co 



wye 



-2/k 



exp 



Co 



ds 



v 2 (c 2 + (3e-^ u )dn 



where u(/i), and in the integrands are obtained from the e = combustion 
wave discussed in Section [2~T1 The apparent dependence of ci on the wave coordinate 
£ is spurious: if I is an anti-derivative of the inner integrals in (|2.16|) . a multiplicative 
term exp [—-?(—£)] emerges in both the numerator and denominator, which therefore 
cancels. Hence, any convenient value for £ can be chosen in (|2.16p . for example 0. 

Equation (|2.16p is a powerful expression in which the wavespeed correction is 
expressed purely in terms of the (unperturbed) infinite Lewis number wavefront and 
system parameters. This correction was obtained through an application of the slow 
manifold and Melnikov's method (suitably modified). While developed within the 
current specific context, we note that this technique can in fact be used in a variety 
of instances which are modeled through coupled reaction-diffusion equations in which 
the diffusivities are very different. 

We note that v < for the infinite Lewis number wavefront, as is clear from 
the phase-portrait Figure 12.11 Alternatively, u is smaller at the front of the wave, 
where fuel is yet to be burnt, and is therefore a decreasing function of /i, leading to 
v = u' < 0. Based on this, (|2 . 1 6[) immediately displays that c\ < 0, proving the 
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property that the wavespeed decreases when fuel diffusivity is included. This is in 
agreement with the numerical observations in Section [2.21 

Equation (|2.16|) provides an explicit perturbative formula on how the wavespeed 
varies through the inclusion of the finitcncss of the Lewis number, expressed entirely 
in terms of the infinite Lewis number combustion wave. This result is used to compute 
the solid line in Figure YTM which is the theoretical wavespeed 0.5707 — 0.1552e ob- 
tained by using (|2.16[) and (|2.12[) when j3 = 1. When e is small, it forms an excellent 
approximation to the numerically obtained wavespeed, as described in Section 12.21 
Indeed, Figure l2~4l show that the theoretical line is tangential to the curve formed by 
the closed circles. 

The perturbation wavespeed C\ as a function of (3 appears in Figure [2~71 There is 
a value of (3 (around 2) at which the absolute influence of the finitcncss of the Lewis 
number is greatest. Nevertheless, since Cq is itself a function of /3, it makes sense to 
investigate the relative influence c\/co of the perturbative term. Such is presented 
in the numerically computed Figure 12.81 The graph is virtually linear and has zero 
intercept. In other words, the complicated quotient in (|2.16p is in fact proportional 
to the unperturbed wavespeed cq{(3), with the proportionality factor independent of 
(3. We therefore arrive at the approximation 

(2.17) c(j3,e) = c (/?) [1 - 0.267 e0\ = c (J3) 

for large Lewis numbers, with excellent validity across all (3, and with Co(/3) also known 
through (|23|) . 

Equation (|2.17p shows that the wavespeed, as a fraction of the infinite Lewis 
number wavespeed, acquires a correction linear in the ratio /3/Le. We are not aware 
of any such result being previously reported in the literature of combustion waves. 
Moreover, the simplicity of this expression is remarkable. For the specific instance 
(3 = 1, we apply this formula in order to arrive at the dashed line in Figure [2Tl Our 
perturbative theory has clearly given us a very accurate and simple approximation, 
and elucidates the straightforward dependence of the wavespeed on the parameters (3 
and Le. 

3. Stability analysis. In this section, we test the stability of the combustion 
wavefront (w, y) = (ito(£), yo(£)) we have found as a solution to (jl.l[) at large Lewis 
numbers. Consider a perturbation of the form 



0.267 — 
Le 



(3.1) 



u = u a (0 + U(£)e xt , y = yo{Q + Y(£)e 
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Fig. 2.8. Relative size of perturbative wavespeed as a function of j3 

At first order, U and Y satisfy an eigenvalue problem 

( U \ 

V 

Y 
V z ) 

Linear instability occurs if there are values of A in the right half plane for which 
(|3.2p possesses a solution uniformly bounded for all £. It turns out that such values of 
A can be investigated by analyzing the Evans function [17], which is a complex analytic 
function E(X) whose zeros correspond to exactly these A values. If for example it can 
be shown that E{X) has no zeros in the right half plane, the indications from the 
point spectrum of (|3.2p is that the wavefront is stable. If there exist zeros of E(X) in 
the right half plane, the wavefront is unstable. A description of the Evans function 
as used in our study is given in Appendix iBl This was proposed in [2, 16,23,40] and 
has been used by [20] for a detailed numerical analysis of (|1.1[) . (It must be also 
mentioned that in the linear stability analysis, it is necessary to consider the essential 
spectrum associated with p.2[) ; it turns out that this has no intersection with the 
right half plane and therefore need not worry us further.) 

We begin with a traveling wave solution uo(£) and t/o(0 obtained using standard 
shooting methods in Section [2j and then compute the Evans function using the proce- 
dure outlined in Appendix [B] We note that the system is very sensitive due to its stiff- 



(3.2) 



/ V \ 

V 
Y 

V z J 



yo e -l/uo 



£yo e -l/«o 
Tuf. 





_ e -l/«o 



A + P e -l/«o 



\ 


1 
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ness. We found a solution to be accurate enough if we obtain E(X = 0) ~ O(10~ 12 ). 
We are guided in our calculations by the detailed stability analysis of Gubernov et 
al [20]. They show, for example, the lack of any eigenvalues of positive real part for 
small (3, but show that for (3 large enough, two eigenvalues pop into the right half 
plane exhibiting a Hopf bifurcation. Physically, this corresponds to a pulsating in- 
stability in the wavefront, a well known phenomenon also occurring for Le = oo even 
in slightly different models [5,14,28,37]. Gubernov et al extend these infinite Lewis 
number analyses by producing in Figure 5 in [20] the stability boundary in (3-e space 
(their r is our e). We verify here that our numerically computed wavefronts display 
the characteristics outlined by them. 

In Figure [3~T1 we show the Evans function E(A) as it varies with increasing A G M 
for Le = 17 and (3=1 (this corresponds to a stable regime in Figure 5 of [20]). We 
find that Evans function does not have any positive real roots. To test for complex 
roots we vary A S iK; using Cauchy's theorem we can calculate the winding number 
to detect possible oscillatory instabilities. In the left panel of Figure 13.21 we show 
the complex Evans function. Since the system (jl.ip is translationally invariant, the 
Evans function has at least a simple zero at A = 0. We checked with a little off-set of 
the order 0(1O -5 ) whether the (real) value of the Evans function at A = is shifting 
towards larger values or smaller values. The off-set allows us to integrate parallel 
to the imaginary axis of A and therefore excluding the zero of the Evans function 
stemming from the root at A = 0. This enables us to attribute roots of the Evans 
function to either the translational mode or to a real instability. For the case depicted 
in Figure I3~2l we find that the Evans function moves to the right. This means that 
the Evans function can be cast in the topologically equivalent form depicted in the 
right panel of Figure [3~2l and it has clearly a winding number zero. We therefore find 
that at these parameter values there are no unstable eigenvalues. (Note that for this 
argument to work we need our definition of the Evans function to be analytic which 
excludes standard methods such as Gram-Schmidt orthogonalizations.) 

We next choose Le = 100 and (3 = 9, parameters at which (according to Figure 5 
in [20]) an oscillatory instability is to be expected. In Figure [3~3l we show the Evans 
function in this situation, with a zoom in displayed in Figure 13.41 To determine the 
winding number we need to check whether the small circle of the Evans function in 
Figure I3~4l includes the zero or not. We can do so by allowing again for a small off- 
set of A. We find that the circle includes zero. Unfolding the behavior of the Evans 
function then allows to sketch a topologically equivalent Evans function as in the right 
panel of Figure (33 We verified that the instability is indeed oscillatory by examining 
the Evans function for A S M + , which reveals no zeros. Hence our wavefront displays 
the predicted characteristics of [20]. Stability properties are not affected unduly by 
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FlG. 3.1. The real part of the Evans function E(X) for Le = 17 and = 1 as a function of A. 
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K(E(\)) 

Fig. 3.2. Left: The real versus the imaginary part of the Evans function E(X) for Le = 17 and 
(3 = 1. The spectral parameter A varies along the imaginary axis. Right: A sketch of a topologically 
equivalent Evans function. The winding number is clearly zero indicating stability. 

the finiteness of the Lewis number, despite the singularity of this limit. 

4. Concluding remarks. In this article, we have studied combustion wavefront 
in a one-dimensional medium. Our concentration was on very high Lewis numbers 
relevant to high-density supercritical combustion. We determine the wavespeed as a 
function of the exothermicity parameter (3 and the Lewis number Le, by seeking the 
wavespeed value which establishes a connection between the fixed points correspond- 
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Fig. 3.3. The real versus the imaginary part of the Evans function E(X) for Le = 100 and 
(3 = 9. The spectral parameter A varies along the imaginary axis. 
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Fig. 3.4. Left: Closeup of the Evans function depicted in Fig. \3.3\ into the region E(X) = 0. 
Right: A sketch of a topologically eguivalent Evans function. The winding number is clearly two 
indicating an oscillatory instability. 



ing to the fully burnt and the unburnt states. The infinite Lewis number instance 
reveals an exponential dependence of the wavespeed on /?, for which we determine 
an empirical formula. We then use several suitably modified dynamical systems tech- 
niques (slow manifold reduction, and Melnikov's method) to compute an explicit 
formula (|2. 16p for the correction to the wavespeed when including the effect of large, 
but not infinite, Lewis number. We hence obtain a simple formula (|2.17p which shows 
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that the relative change in the wavespeed is proportional to /3/Le for large Lewis 
numbers. Our theory is shown to have excellent consistency with the numerically 
computed wavespeed for large Le, as we show in Figure \TM 

The stability of the high Lewis number wavefronts is then tested numerically 
based on the Evans function technique. Our results are in agreement with the stability 
boundaries presented by Gubernov et al in Figure 5 in [20] . 

We remark that the modified Melnikov's method that we have used can in fact be 
used in more general situations which are described by two coupled reaction-diffusion 
equations with strongly differing diffusivities. Based on a known wavefront or wave- 
pulse solution for when the smaller of the diffusivities is zero, our technique can in some 
instances be used to determine the wavespeed correction resulting from the inclusion 
of the (previously neglected) diffusivity. Alternatively, it can be adapted to situations 
in which the wavespeed changes due to some other small parameter. Our analysis of 
high Lewis number wavefronts therefore provides a new perturbative methodology for 
analyzing certain classes of reaction-diffusion equations, pattern formation problems 
and combustion waves. 
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Appendix A. Melnikov function derivation. 

We briefly outline the modifications needed to the standard Melnikov approaches 
[4,21,38] relevant to Section \2A\ Our system is z = / (z) + eg (z), as given in (|2.13|) . 
Consider a particular parametrization of the heteroclinic z(£). Imagine the perturbed 
system as embedded in three-dimensional (z, s) space. In a "time"-slice s = so, let 
T be the normal vector to the heteroclinic drawn at the point z(0) = zq. The usual 
approach is to compute the distance between the perturbed manifolds measured along 
T, and this is expandable as 

(A.I) d( S0 ,e) = e^ + O(e*). 

I f ( z o)| 

Let z u (s) be the trajectory of the perturbed flow which intersects T and which back- 
wards asymptotes to the perturbed fixed point a(e). In other words, z u (s) is a trajec- 
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tory lying on a(e)'s unstable manifold. The standard approach [4,21] is to represent 

z CT (s) =z(s-s )+ez^(s) + C>(e 2 ) 

where a = u (for "unstable"), valid for — oo < s < sq . A similar expansion on so < 
s < oo with a — s (for "stable"), works for the trajectory z s (s), which intersects T 
on the time-slice so and which lies on the stable manifold of the perturbed fixed point 
b(e). Then, the standard Melnikov development (see [4,21]) allows the representation 

A«(s )-A s (s ) 



(A.2) 
where 



d(s ,e) 



If Ml 



<D{e l 



A'( S ) = f(z( S - S „))AzJ( S ) 
for a = u and a — s. Now, [4,21] derive that 

(A.3) A CT = V • f (z(s - s )) A" + f - s )) A g (z(s - s ),s) + 0(e) 

but since the unperturbed dynamical system is volume-preserving, have the luxury 
of ignoring the first term on the right hand side. We cannot do so here, but we can 
neglect the second argument in g since our case is autonomous. To deal with the first 
term, we multiply (|A.3j) by the integrating factor 



/i(s) = exp 



[ V-f (z(r-s )) dr 
Jo 



before proceeding. Having done so, we integrate from — oo to sq by choosing a = u, 
then integrate from sq to oo by choosing a = s, and then add the two equations to 
get 



a«( S0 )-a s ( S0 ) = r 

J —< 



M(so) 



f Ag(z(s - s )) ds. 



(This is an adaptation of the standard process [4,21].) In conjunction with (|A.1[) 
and (|A.2j) . and also employing the shift s — s — > s in the integrand, we obtain the 
Melnikov function 



M (so) = 



/oo r ps 

exp — / V-f (z(r)) dr 
-oo L Jo 



f A g (z(s)) ds 



which no longer depends on so. Having dealt with the non volume-preserving instance, 
the next step is to change our attitude: rather than measuring the distance in a time- 
slice s but at a specific point Zq, we ignore time-slices (since our perturbed system 
is itself autonomous), and allow the point to vary along the heteroclinic. To do so, 
choose a different parametrization w(s) = z(s — £) of the heteroclinic. Thus, the point 
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w = w(0) = z(— £) can be varied along the heteroclinic by choosing different values 
of £. Therefore, £ will represent different points along the heteroclinic at which the 
distance measurement is to be made (cf. "heteroclinic coordinates" of [38]). Using 
the w trajectory, our earlier results can be expressed as 



(A.4) 
where 



|f(wo)| 



|fW-0)l 



0(e'), 



M(f) 



exp 



exp 



/ V • f (w(r))dr 
Jo 

V • f (z(r - £)) dr 



f A g (w(s)) ds 

f Ag(z(s-£))ds 

f A g (z(s)) (is 



exp - / V • f (z(r - £)) dr 

- oo JO 

/oo r /* s 

exp -/ V-f(z(r))dr fAg(z(s))ds. 
-oo L ./-£ 

This, in conjunction with (|A.4|1 . is the expression used in Section 
Appendix B. Evans function definition. 

Here, we describe the Evans function approach for analyzing linear stability. In 
general, the linear stability of a localized traveling wave solution to a system of PDEs 
is obtained by studying the eigenvalue problem 



(B.l) 



Cw — Xw , 



where the matrix differential operator C arises from the linearization of the PDEs. 
The traveling solution is said to be linearly stable if the spectrum of C lies in the 
closed left half-plane. 

The system (|B.1[) can be turned into a linear dynamical system of the form 



(B.2) 



U t = A(£, A) U 



where A(£, A) is an n x n square matrix depending on £ = x — ct and the spectral 
parameter A (in our case, n = 4). It can be shown that the asymptotic behavior of 
the solutions to (|B.2[) is determined by the matrices 



L±oo 



(A) 



lim A(£, A) 

J— »±oo 



in the following sense (see [11] for details). If /i + (resp. is an eigenvalue of A +00 
(resp. A-_oo) with eigenvector v + (resp. v~), then there exists a solution w + (resp. 
w~) to (|B.2[) with the property that 



(B.3) 



lim w^e M ^ 
6 — 



resp. lim w e 

£ — oo 
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Note that the superscript "+" refers to exponentially decaying behavior at £ = +00, 
while "— " refers to £ = —00. 

To study the linear stability, one should consider both the essential and point 
spectrum of C. The essential spectrum of C consists of the values of A for which Aoo 
or A-oo has purely imaginary eigenvalues [22]. The point spectrum can be studied by 
means of the Evans function, first introduced by Evans [16] and later generalized [2]. 
Roughly speaking, the zeros of this complex-valued function are arranged to coincide 
with the point spectrum of C. 

Let fl denote a domain of the complex A plane with no intersection with the 
essential spectrum and let n s and n u denote, respectively, the number of eigenvalues 
of Aqo with negative real part and the number of eigenvalues of A-oo with positive real 
part in fi. We assume that n s +n u = n. Let wf(\,£), i = l,2,...,n s (resp. w~(A,£), 
i = 1,2,..., n u ) be linearly independent solutions to (1B.2|) converging to zero as £ — * 00 
(resp. £ — » —00) which are analytic of A in Q. Clearly, a particular value of A belongs 
to the point spectrum of £ if (|B.2[) admits a solution that is converging to zero for 
both £ — + ±00, that is if the space of solutions generated by the wf intersects with the 
one generated by the w~ . To detect such values of A in fi, one can use the definition 
of the Evans function given in [33] 



in which the w~ are evaluated at £ = 0. This function is analytic in fi, is real for real 
values of A and the locations of the zeros of E(X) correspond to eigenvalues of C 

The first numerical computation of the Evans function was by Evans himself 
in [17], and followed by [32,35]. However, in all three papers n s = 1, in which case a 
standard shooting argument can be used. In standard shooting algorithms one follows 
the stable and/or unstable manifolds at £ — ±00. The Evans function is then given 
as the intersection of these manifolds. As shown in Section [3j our system has n — 4 
and n s = n u — 2. This causes the following practical problem: although the n s (or n u 
respectively) eigenvectors are linear independent solutions of the eigenvalue problem 
(|B.2j) at £ = ±00, the numerical integration scheme will lead to an inevitable alignment 
with the eigen-direction corresponding to the largest eigenvalue. This collapse of 
the eigen-directions is usually overcome by using Gram-Schmidt orthogonalization. 
However, this is not desirable for calculating the Evans function as it is a nonanalytic 
procedure which then subsequently prohibits the use of Cauchy's theorem (argument 
principle) to locate complex zeros of the Evans function. The Evans function is 
therefore best calculated using exterior algebra [1,3,9,11,12,15,31,34]. 

We briefly review the method here, with specific regard to the situation in which 
n = 4 and n s = n u = 2. For more details the reader is referred to [1,3,11,15], and to 
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the numerical computation in [20]. The main idea behind exterior algebra methods 
(or compound matrices methods) is that the linear system (|B.2|) induces a dynamical 
system on the wedge-space /\ 2 (C 4 ) for n s = n u = 2. The wedge-space /\ 2 (<C 4 ) is the 
space of all two forms on C n . This is a space of dimension ( 4 ) =6. The induced 
dynamics on the wedge-space /\ 2 (C 4 ) can be written as 

(B.4) U^A( 2 )(£)U, UeA 2 (C 4 ). 

Here the linear operator (matrix) A^ 2 ) is the restriction of A(£, A) = {a^ } to the 
wedge space /\ 2 (C 4 ). Using the standard basis of /\ 2 (C 4 ) 

,„ w 1 =e 1 Ae 2 , w 2 = eiAe 3 , uj 3 =e 1 /\e 4 , 

(B.5) 

W4 = e 2 Ae 3 , w 5 = e 2 Ae 4 , W6 = e 3 Ae 4 , 

where ex 2,3,4 is the standard basis of C™, we can find the matrix A^ 2 ) : A 2 (C 4 ) 
/\ 2 (C 4 ) as a complex 6x6 matrix. With respect to the basis (|B.5[) . A^ 2 ) takes the 
explicit form 



0.11 +0,22 


a.23 


a24 


— O13 


— ai4 





032 


«ll+ffl33 


034 


a.12 





-O14 


a.42 


O43 


Oil +O44 





012 


Ol3 


—031 


031 





022+033 


034 


-024 


— a« 
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a.43 


022+044 


023 





— a 4 l 


031 


— 042 


032 


033+044 



General aspects of the numerical implementation of this theory and details for these 
constructions in more general systems can be found in [3, 11]. 

Linearity assures that the induced matrix A( 2 )(£,A) is also differentiable and 
analytic. Hence, the limiting matrices, 

Ail(A)= Urn A< 2 >(£,A), 
£-+±00 

(2) 

will exist. It can readily be shown that the eigenvalues of the matrix A^^A) con- 
sists of all possible sums of 2 eigenvalues of A± OQ (A). Therefore, for 5ft(A) > 0, the 

eigenvalue of A^ | 2 ] (A) with the most negative real part is given by o+(A) = /i^ + [i\ . 

( 2) 

The eigenvalue o+(A) has real part strictly less than any other eigenvalue of A y + ^ X1 (X). 

Analogously, the eigenvalue of A_4 ) (A) with the largest non- negative real part is given 

by cr-(A) = + ^2 ■ The eigenvalue c_(A) has real part strictly greater than any 
(2) 

other eigenvalue of A_4>(A). Note that the eigenvalues a± are simple, and analytic 
functions of A. 

Let £ (A) be the eigenvectors associated with o±(A), defined by 
(B.6) aV1(A)C + (A) = a + (AR + (A) and A^l(X)C(X) = <r-(A)C(A) . 
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These vectors can always be constructed in an analytic way (see [11]), and are readily 
found to be C ± (A) = vf A vf. 

Let U ± (^,A) G /\ 2 (C 4 ) be the solution of the linear system (|B.4j) which satisfy 
lim^^ioo e _cr± ( A ^U ± (^, A) = C ± (^)- This allows us to define the Evans function as 

(B.7) £(A) =7VU-(e,A) AU+K,A), A e A , 

where 

(B.8) N = e~ & r(s,\)ds and r(^, A) = Tr(A(£, A)) . 

Expressing U ± (^, A) as a linear combination with respect to the basis (IB.5|) 

6 
3 

the expression for the Evans function (|B.8|) can be simplified to 
(B.9) £(A)=AqU-(£,A),£U+(£,A)] 6 , 

where [•, -]6 is the complex inner product in C , and the representation of the Hodge- 
star operator £ in the basis (|B.5|) is 

1 " 
0-10 

1 





Using the Hodge-star operator, we can relate the the most unstable solution U~ of the 
linearized system at £ = — oo with the most unstable solution of the adjoint system 
of (|B.4[) at £ — — oo. Details can be found in [3,10,11]. This suggests a normalization 
of the asymptotic eigenvectors according to 

(B.io) [C,£C+] 6 = i, 

which assures that E(X) — ► 1 for |A| — ► oo. 

Note that the translational invariance of (|1.1|) guarantees that the Evans function 
can be evaluated at any (fixed) spatial location £*. However, to avoid unwanted 
growing of the solutions U ± we will consider the scaled solutions 



£ = 



























Q 








1 





-1 





1 









(B.ll) 



1^(6 A) = e- <T±(A )«U ± (£,A) 
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The scaling (|B.11[) ensures that U + (£, A)|^_, 4 is bounded. The corresponding equa- 
tion on A 2 (C 4 ) 

(B.12) = [AW(£,X)-a ± (\)I]lJ±, U±(£,A)| s=i±oo =C ± (A), 

is integrated from £ = L± oc to £ = £* (where is arbitrary but fixed). 

The system (|BT2|) can be integrated using the second-order implicit midpoint 
method. For a system in the form = B(£, A)U, each step of the implicit midpoint 
rule takes the form 

(B.13) U" +1 = [I - + ±AxB„ +1/2 ] U" , 

where B n+1/2 = B(x n+1/2 , A). 
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